/*  partition.cpp
 * Copyright (C) 2010, Francisco Claude & Patrick Nicholson, all rights reserved.
 *
 * Francisco Claude & Patrick Nicholson <(fclaude|p3nichol)@cs.uwaterloo.ca>
 *
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <stdio.h>
#include <string>
#include <set>
#include <iostream>
#include <vector>
#include <map>
#include <list>
#include <limits>
#include <cmath>
#include <cstdlib>
#include <fstream>

using namespace std;

/**
 * A simple class for storing 2D points.  Complete with overloaded
 * comparator and stream output.
 */
class Point2
{
  public:

    double x, y;
    Point2(const Point2 &p) { x = p.x; y = p.y; }
    Point2(double _x, double _y) { x = _x; y = _y; }
    Point2() { x = 0; y = 0; }

    bool operator <(const Point2& b) const
    {
      return x < b.x || (x == b.x && y < b.y);
    }

    friend ostream& operator<< (ostream &out, Point2 &cPoint);
};

ostream& operator<< (ostream &out, Point2 &p) {
  out << "(" << p.x << "," << p.y << ")";
  return out;
}


/**
 * Layer class, used to encapsulate the creation of layer structures
 * on the input data set, and also allow the extraction of the longest
 * increasing subsequence from the data.  Currently only decending
 * layers are supported.
 */
class Layers
{

  private:

    int n;
    vector<list<Point2> > layerlist;

  public:

    /**
     * The Layer constructor takes an STL set of points as input, and
     * then builds the set of layers using the algorithm described in
     * Yehuda and Fogel's paper.
     * Assumes thePoints is non-empty.
     */
    Layers(set<Point2> *thePoints) {

      Point2 *P = new Point2[thePoints->size()];
      n = 0;

      for (set<Point2>::iterator sit = thePoints->begin(); sit != thePoints->end(); sit++) {
        P[n++] = *sit;
      }

      double *A = new double[n+1];
      A[0] = numeric_limits<double>::max();

      int l = 0;
      for (int i = n-1; i >= 0; i--) {
        int j = 0;
        int j_prime = 0;
        for(; j <= l; j++) {
          if (A[j] > P[i].y) {
            j_prime = j;
          }
        }
        j = j_prime;

        if (j == l) {
          layerlist.push_back(list<Point2>());
          l++;
        }

        list<Point2> *L = &layerlist[j];
        L->push_front(Point2(P[i]));
        A[j+1] = P[i].y;
      }
      delete [] P;
      delete [] A;
    }

    ~Layers() {
      layerlist.empty();
    }

    int size() {
      return layerlist.size();
    }

    list<Point2> getLongestIncreasingSubsequence(int k) {

      list<Point2> lis;

      if (layerlist.size() > 0) {
        int l = layerlist.size() - 1;
        Point2 q_prime = layerlist[l].front();
        lis.push_front(Point2(q_prime));
        for (int i = l - 1; i >= l - k; i--) {
          for (list<Point2>::iterator lit = layerlist[i].begin(); lit != layerlist[i].end(); lit++) {
            Point2 q = *lit;
            if (q.x >= q_prime.x && q.y >= q_prime.y) {
              lis.push_front(Point2(q));
              q_prime = Point2(q);
              break;
            }
          }
        }
      }
      return lis;
    }
};

bool descending;

/**
 * Partition the set "points" into two sets, "redPoints" and
 * "bluePoints", using (a not-so-efficient-version of) the algorithm
 * described in Yang et al. 2007 "A Comparative Sutdy of Efficient
 * Algorithms for Partitioning a Sequence into Monotone Subsequences."
 * It is less efficient since instead of using the layer update method
 * described by Yehuda and Fogel, it just recomputed the layers at
 * each step in the algorithm.  Thus, the complexity is
 * O(n*lg(n)*sqrt(n)), instead of O(n*sqrt(n)).
 */
void computePartition(set<Point2> *points, set<Point2> *redPoints, set<Point2> *bluePoints) {

  if (points->empty()) {
    return;
  }

  redPoints->clear();
  bluePoints->clear();

  for (set<Point2>::iterator sit = points->begin(); sit != points->end(); sit++) {
    redPoints->insert(*sit);
  }

  int sequences = 0;
  int K = (int) (sqrt(2.0 * points->size() + 0.25) - 0.5);

  for (int k = K; k > 0; k--) {
    Layers redLayers = Layers(redPoints);

    int l = redLayers.size();
    //cout << "k=" << k << " l=" << l << endl;
    if (l <= k) {
      sequences += l;
      break;
    }
    else {
      list<Point2> lis = redLayers.getLongestIncreasingSubsequence(k);
      for (list<Point2>::iterator it = lis.begin(); it != lis.end(); it++) {
        Point2 p = *it;
        bluePoints->insert(p);
        redPoints->erase(p);
      }
      sequences++;
      if (redPoints->size() == 0) {
        break;
      }
    }
  }

  cerr << "Number of points:    " << points->size() << endl;
  cerr << "Number of sequences: " << sequences << endl;
  cerr << "Bound on sequences:  " << K << endl;
}


/**
 * Read the input file, which contains lines of the format:
 *
 * <point number> <x-coordinate> <y-coordinate>
 *
 * where the x and y coordinates are floating point numbers, and dump
 * the points into an STL set. Returns true if the file as read
 * without errors.
 *
 */
bool loadPoints(string inFilePath, set<Point2> *points) {
  ifstream infile (inFilePath.c_str());
  string line;
  int lineno = 1;
  if (infile.is_open()) {
    while (! infile.eof() ) {
      getline (infile,line);
      double px, py;
      int match = sscanf(line.c_str(), "%*d %lf %lf", &px, &py);
      if (match == 2) {
        if(descending)
          points->insert(Point2(px,py));
        else
          points->insert(Point2(-px,py));
      }
      else {
        //cerr << "Error parsing file on line " << lineno << ": " << line << endl;
        //cerr << "Only matched " << match << " tokens. Skipping line." << endl;
      }
      lineno++;
    }
    infile.close();
    return true;
  }
  return false;
}


/**
 * Write the points from redPoints and bluePoints into a file using
 * the same format as the input file, except adding a second column:
 *
 * <point number> <ascending/descending> <x-coordinate> <y-coordinate>
 *
 * Returns true if the file is created and written to without errors.
 */
bool writePartitionedPoints(string outFilePath, set<Point2> *redPoints, set<Point2> *bluePoints) {

  ofstream outfile(outFilePath.c_str());
	outfile.precision(14);
  int pNum = 0;
  if (outfile.is_open()) {

    if(descending) {
      outfile << redPoints->size() << " " <<  bluePoints->size() << endl;
      for (set<Point2>::iterator it = redPoints->begin(); it != redPoints->end(); it++) {
        Point2 p = *it;
        outfile << pNum++ << " D " << p.x << " " << p.y << endl;
      }

      for (set<Point2>::iterator it = bluePoints->begin(); it != bluePoints->end(); it++) {
        Point2 p = *it;
        outfile << pNum++ << " A " << p.x << " " << p.y << endl;
      }
    }
    else {
      outfile << bluePoints->size() << " " <<  redPoints->size() << endl;
      for (set<Point2>::iterator it = bluePoints->begin(); it != bluePoints->end(); it++) {
        Point2 p = *it;
        outfile << pNum++ << " D " << -p.x << " " << p.y << endl;
      }

      for (set<Point2>::iterator it = redPoints->begin(); it != redPoints->end(); it++) {
        Point2 p = *it;
        outfile << pNum++ << " A " << -p.x << " " << p.y << endl;
      }
    }
    outfile.close();
    return true;
  }
  else {
    return false;
  }
}


/**
 * Main function: read the input file of points, partition them using
 * the computePartition function, then write the output to the
 * specified output file.
 */
int main(int argc, char **argv) {

  set<Point2> points, redPoints, bluePoints;

  if (argc != 4) {
    cerr << "Usage: " << argv[0] << " <infile> <outfile> <direction>" << endl;
    exit(1);
  }

  descending = (string(argv[3])==string("descending"));
  cout << "Direction: " << descending << endl;
  string inFilePath(argv[1]);
  string outFilePath(argv[2]);

  if (loadPoints(inFilePath, &points)) {
    cerr << "Opened file \"" << inFilePath << "\"" << endl;
    computePartition(&points, &redPoints, &bluePoints);
    if (writePartitionedPoints(outFilePath, &redPoints, &bluePoints)) {
      cerr << "Wrote output to file \"" << outFilePath << "\"" << endl;
    }
    else {
      cerr << "Failed to write output to file \"" << outFilePath  << "\"" << endl;
      return 1;
    }
  }
  else {
    cerr << "Error opening file \"" << inFilePath << "\" for reading." << endl;
    return 1;
  }
  return 0;
}

